 ***THIS PROGRAM EVALUATES FGLS (PARKS);
libname PARKSSUR 'E:\MANTO\UC\THESIS\CHAPTER 3\ARTICLE\REPLICATION FORLER\PARKSREP';

proc iml;
start program;
use PARKSSUR.tourismitaly;
read all var {lcrime} into y1;
read all var {growth, lgdp, lun, ldens, ltourism, ldiplomapc, ldetcrime} into x1;


T = 18;
n = 17;
nT = n*T;
reps = 500;
nboots = 499;
M0 = 1;
y = j(nT,1,0);
y = y1[1:nT];
Ones = j(95*T,1,1);
x2 = Ones||x1;
X= j(nT,ncol(x2)*n,0);
R = j(1,ncol(x2)*n,0);
R[2] = 1;

R2 = j(1,ncol(x2)*n,0);
R2[2] = 1;
R2[10] = -1;


R3 = j(2,ncol(x2)*n,0);
R3[1,2] = 1;
R3[1,10] = -1;
R3[2,3] = 1;
R3[2,11] = -1;

seed = 4060;

** IMPLEMENT SUR;
X = j(nT,ncol(x2)*n,0);
do i = 1 to n;
        X[(1+(i-1)*T):i*T,(ncol(x2)*(i-1)+1):ncol(x2)*i] = x2[(M0-1)*T+(1+(i-1)*T):(M0-1)*T+i*T,]; 
end;

start ParksSur(rhohatsur, Test, sigmahat, H, A, P, bparks, vcovbparks, T, y, X); 
        n = nrow(X)/T;
        nT = n*T;
        
        ** OLS Residuals;

        eols = y-X*solve(X`*X,X`*y);
        ee = j(T,n,0);
        do i = 1 to n;
                ee[,i] = eols[1+(i-1)*T:i*T];
        end;
        sigmahatsur = (1/T)*(ee`*ee);

        start MyInv1(InvM, M);  
                D0 = sqrt(diag(M));
                R0 = ginv(D0)*M*ginv(D0);
                InvM = ginv(D0)*ginv(R0)*ginv(D0);
        finish MyInv1;

        run MyInv1(invsigmahatsur, sigmahatsur);                
        invomegahatsur = invsigmahatsur@I(T);   

        start MyInv2(InvM, M);
                call SVD(Ustar, qstar, Vstar, M);
                InvM = Vstar*ginv(diag(qstar))*Ustar;
        finish MyInv2;
        
        v = X`*invomegahatsur*X;
        run MyInv1(invv, v);
        bsur = invv*(X`*invomegahatsur*y);
        esur = y - X*bsur;
        
        *** COMPUTE PARKS ESTIMATOR USING THE SUR RESIDUALS;
        eesur = j(T,n,0);
        do i = 1 to n;
                eesur[,i] = esur[1+(i-1)*T:i*T];
        end;
        rhohatsur = j(n,n,0);
        rhohati = j(n,1,0);
        do i = 1 to n;
                eesur1 = eesur[2:T,i];
                eesur2 = eesur[1:T-1,i];
                rhohatsur[i,i] = eesur1`*eesur2/(eesur2`*eesur2);               
        end;

        *** find P0;

        P0 = j(nT,nT,0);
        ** I leave zeros at for the first observation. I will cut those below and upgrade this matrix later to P;
        do i  = 1 to n;
                do j = 2 to T;
                        P0[(i-1)*T + j,(i-1)*T + j] = 1;
                        P0[(i-1)*T + j,(i-1)*T + j-1] = -rhohatsur[i,i];
                end;
        end;

        *** Transform y and x using P0;
        y_trans = P0*y;
        x_trans = P0*X;
        yreduced = j(nT-n,1,0);
        xreduced = j(nT-n,ncol(X),0);
        do i = 1 to n;
                yreduced[1+(i-1)*(T-1):i*(T-1)] = y_trans[2+(i-1)*T:i*T];
        xreduced[1+(i-1)*(T-1):i*(T-1),] = x_trans[2+(i-1)*T:i*T,];
        end;
        /*  reduce transfoormed Y and X by deleting the first values for
     each equation then run an OLS model and compute sigmahat       */;
        breduced = solve(xreduced`*xreduced, xreduced`*yreduced);
        ereduced = yreduced - xreduced*breduced;
        eereduced = j(T-1,n,0);
        do i = 1 to n;
                eereduced[,i] = ereduced[1+(i-1)*(T-1):i*(T-1)];
        end;
        sigmahat = (1/(T-1))*eereduced`*eereduced;

        *** compute V0;
        V0 = j(n,n,0);
        do i = 1 to n;
                do j = 1 to n;
                        V0[i,j] = sigmahat[i,j]/(1 - rhohatsur[i,i]*rhohatsur[j,j]);
                end;
        end;
        /*
        eigenval=eigval(V0)||eigval(sigmahat);
        
        if ncol(eigenval) > 2 then Test  = 0;
        if ncol(eigenval) = 2 then do;
                m = min(eigenval);       
                if m <= 0 then Test = 0;
                if m > 0 then Test = 1;
        end;
        
        
        if Test = 0 then goto TheEnd;
        */

        **** Compute A;  
        
     H = root(sigmahat, "NoError");
        B = root(V0, "NoError");
                Test = 1;
        if any(B = .) then Test = 0;
        if any(H = .) then Test = 0;
        if Test > 0 then do;
 			detB = det(B/B[:]);
			if detB = 0 then Test = 0;
			if Test = 1 then do;
                *run MyInv2(Bi, B`);
				Bi = inv(B`);
                A = H`*Bi;
                *** Construct P by completing P0;
                P = P0;
                do i = 1 to n;
                        do k = 1 to i;
                                P[1+(i-1)*T, 1+(k-1)*T] = A[i,k];
                        end;
                end;

                *** Use P to transform Y and X;
                Ystar = P*y;
                Xstar = P*X;

                *** Use transformed Y and X to compute Parks estimator and coefficient standard errors;         
                run MyInv1(invsigmahat, sigmahat);
                v1 = Xstar`*(invsigmahat@i(T))*Xstar;
                run MyInv1(invv1, v1);  
                bparks = invv1*(Xstar`*(invsigmahat@i(T))*Ystar);
                vcovbparks = invv1;
            end; 
		end;
finish ParksSur;


start PCSE(vstar, bstar, covb, z, y, n, t);

start MyInv1(InvM, M);  
          D0 = sqrt(diag(M));
          R0 = ginv(D0)*M*ginv(D0);
          InvM = ginv(D0)*ginv(R0)*ginv(D0);
finish MyInv1;

nt = n*t;
run MyInv1(invz, z`*z);
b=invz*z`*y;
e=y-z*b;

*****USE OLS RESIDUALS TO CALCULATE COMMON AUTOCORRELATION PARAMETER*********;
ee=j(t,n,0);
do i=1 to n; 
   ee[,i]=e[((i-1)*t)+1:i*t];
end;
rhohati=0;
do i=1 to n;
   rhatnum=ee[2:t,i]`*ee[1:(t-1),i]; 
   rhatden=ee[1:(t-1),i]`*ee[1:(t-1),i];
   check=(rhatnum/rhatden);
   rhohati=rhohati+check;
end;
rhohat=rhohati/n;

if (rhohat < 1) & (rhohat > -1) then do;

***CONSTRUCT THE TRANSFORMATION MATRIX PSTAR*****;
pstar=j(nt,nt,0);
pstari=j(t,t,0); 
pstari[1,1]=sqrt(1-(rhohat**2));
do i=2 to t;
   pstari[i,(i-1)]=-rhohat;
   pstari[i,i]=1;
end;
do i=1 to n;
   pstar[((i-1)*t)+1:i*t,((i-1)*t)+1:i*t]=pstari;
end;

***USE RESIDUALS FROM TRANSFORMED EQUATION TO CALCULATE CROSS-SECTIONAL COVARIANCES**;
zstar=pstar*z;
ystar=pstar*y;

/* 
start MyInv1(InvM, M);  
          D0 = sqrt(diag(M));
          R0 = ginv(D0)*M*ginv(D0);
          InvM = ginv(D0)*ginv(R0)*ginv(D0);
finish MyInv1; 
*/

run MyInv1(vinv, zstar`*zstar);
bstar = vinv*zstar`*ystar;
estar=ystar-(zstar*bstar);
eestar=j(t,n,0);
do i=1 to n; 
   eestar[,i]=estar[((i-1)*t)+1:i*t];
end;
sigma=(1/t)*eestar`*eestar;
vstar=sigma@i(t);
****CALCULATE PANEL-CORRECTED STANDARD ERRORS*****;
run MyInv1(invzstar, zstar`*zstar);
covb=invzstar*zstar`*vstar*zstar*invzstar;
seb=sqrt(vecdiag(covb));
end;

finish PCSE;

run ParksSur(rhohatsur, Test, sigmahat, H, A, P, bparks, vcovbparks, T, y, X);

g1 = (R*bparks)`*inv(R*vcovbparks*R`)*(R*bparks);
g2 = (R2*bparks)`*inv(R2*vcovbparks*R2`)*(R2*bparks);
g3 = (R3*bparks)`*inv(R3*vcovbparks*R3`)*(R3*bparks);

*run PCSE(vstar, bstar, covb, X, y, n, t);


bparksres = bparks;
bparksres[2] = 0;

start delcol(X,i);
	return(x[,setdif(1:ncol(x),i)]);
finish delcol;
Xres = delcol(X,2);
* Xres = X[,1]||X[,3:ncol(X)];

bparksres2 = bparks;
bparksres2[2] = bparks[10];
Xres2 = delcol(X,10);
Xres2[,2] = X[,2] + X[,10];



bparksres3 = bparks;
bparksres3[2] = bparks[10];
bparksres3[3] = bparks[11];
Xres3 = delcol(X,{10 11});
Xres3[,2] = X[,2] + X[,10];
Xres3[,3] = X[,3] + X[,11];

*** THE BOOTSTRAP STARTS HERE;
Testres = 0;
do until (Testres = 1);
	Testrank = 0;
	do until (Testrank = 3);
        v = normal(j(n,T,seed));
        e = H`*v;
        u = j(n, T, 0);
        run MyInv2(InvA, A);
        u[,1] = InvA*e[,1];
        do i = 2 to T;
                u[,i] = rhohatsur*u[,i-1] + e[,i];
        end;
        ysim1 = X*bparksres + shape(u,nT, 1);
        ysim2 = X*bparksres2 + shape(u,nT,1);
        ysim3 = X*bparksres3 + shape(u,nT,1);


        asympt_crit_value1 = 3.841455338;
        asympt_crit_value2 = 3.841455338;
        asympt_crit_value3 = 5.991464547;

        run ParksSur(rhohatr, Testr, sigmahatr, Hr, Ar, Pr, bparksr, vcovbparksr, T, ysim1, Xres);
                *A_boot = Ar;
                *H_boot = Hr;
                *rhohat_boot = rhohatr; 
        run ParksSur(rhohatr2, Testr2, sigmahatr2, Hr2, Ar2, Pr2, bparksr2, vcovbparksr2, T, ysim2, Xres2);
        run ParksSur(rhohatr3, Testr3, sigmahatr3, Hr3, Ar3, Pr3, bparksr3, vcovbparksr3, T, ysim3, Xres3);
		rank1 = round(trace(ginv(Ar)*Ar));
		rank2 = round(trace(ginv(Ar2)*Ar2));
		rank3 = round(trace(ginv(Ar3)*Ar3));
		if rank1 = ncol(Ar) then Testrank = Testrank+1;
		if rank2 = ncol(Ar2) then Testrank = Testrank+1;
		if rank3 = ncol(Ar3) then Testrank = Testrank+1;
	end; 
        Testres = min(Testr, Testr2, Testr3); 
end;

start Statistic(Teststat, bootcrit, bootcrit_np, n, nboots, R, T, ysim, X, bparksres, asympt_crit_value, 
				Xres, Ar, Hr, rhohatr, bparksr, seed); 

	nT = n*T;
	gboot = j(nboots,1,0);
	gboot_np = j(nboots,1,0);


       *transforming the correlated residuals to independent residuals;

        e_original = ysim - Xres*bparksr;
        ee_original = shape(e_original, n, T);
        v_original = j(n,T,0);
        v_original[,1] = Ar*ee_original[,1];
        do i = 2 to T;
             v_original[,i] = ee_original[,i] - rhohatr*ee_original[,i-1];
        end;


       InvHr = ginv(Hr);
       *run MyInv2(InvHr, Hr);
       u_original =InvHr*v_original;
		
       *centering and whitening the inverted residuals;


	u_centered = u_original - u_original[,:];
		
	vcov_u = u_centered*u_centered`;

	*UpperT = root(vcov_u);
	UpperT = root(vcov_u, "NoError");

	UpperT = root(vcov_u, "NoError");
if all(UpperT = .) then Teststat = 0;
else Teststat = 1;
if Teststat = 1 then do;

	u_whitened = ginv(UpperT)`*u_centered;



	boot = 0;
	do until (boot = nboots);
		boot = boot +1;
		*** parametric bootstrap;	
		vp_boot = normal(j(n,T,seed));
		ep_boot = Hr`*vp_boot;
		up_boot = j(n,T,0);

		InvAr = ginv(Ar);
        *InvAr = inv(Ar);

		up_boot[,1] = InvAr*ep_boot[,1];

		do i = 2 to T;
			up_boot[,i] = rhohatr*up_boot[,i-1] + ep_boot[,i];
		end;
		yboot = Xres*bparksr + shape(up_boot,nT,1);
		run ParksSur(rhohatbu, Testboot, sigmahatbu, Hbu, Abu, Pbu, bparksbu, vcovbparksbu, T, yboot, X);
		if Testboot = 1 then do;
			gbu = (R*bparksbu)`*inv(R*vcovbparksbu*R`)*(R*bparksbu);
			gboot[boot] = gbu;
		end;
		
		*** non parametric bootstrap;		
		u_sample = j(n,T,0);
		do i = 1 to T;
			i0 = ceil(T*UNIFORM(i*seed));
			u_sample[,i] = u_whitened[,i0];
		end;


		enp_boot = Hr`*u_sample;


		unp_boot = j(n,T,0);	                
                unp_boot[,1] = InvAr*enp_boot[,1];
                do i = 2 to T;
                        unp_boot[,i] = rhohatr*unp_boot[,i-1] + enp_boot[,i];
                end;
		yboot_np = Xres*bparksr+ shape(unp_boot, nT,1);
		run ParksSur(rhohatbu_np, Testboot_np, sigmahatbu_np, Hbu_np, Abu_np, Pbu_np, bparksbu_np, 
					vcovbparksbu_np, T, yboot_np, X);
		if Testboot_np = 1 then do;
			gbu_np = (R*bparksbu_np)`*inv(R*vcovbparksbu_np*R`)*(R*bparksbu_np);
			gboot_np[boot] = gbu_np;
		end;
		if Testboot = 0 then boot = boot-1;
		if Testboot_np = 0 then boot = boot-1;

		if boot = 0 then boot = boot+1;

	end;
	percentile = .95;
	call qntl(bootcrit, gboot, percentile, 3);
	call qntl(bootcrit_np, gboot_np, percentile, 3);

end;
	
finish Statistic;

*run Statistic(Teststat, bootcrit1, bootcrit_np1, n, nboots, R, T, ysim1, X, bparksres, asympt_crit_value1, Xres, 
                                Ar, Hr, rhohatr, bparksr, seed);
*run Statistic(Teststat, bootcrit2, bootcrit_np2, n, nboots, R2, T, ysim2, X, bparksres2, asympt_crit_value2, Xres2,
                                Ar2, Hr2, rhohatr2, bparksr2, seed);
run Statistic(Teststat, bootcrit3, bootcrit_np3, n, nboots, R3, T, ysim3, X, bparksres3, asympt_crit_value3, Xres3, 
                                Ar3, Hr3, rhohatr3, bparksr3, seed);
/*Stat_crit_values = j(3,4,0);
Statistics = g1||g2||g3;
asympt_crit = asympt_crit_value1||asympt_crit_value2||asympt_crit_value3;
boot_crit = bootcrit1||bootcrit2||bootcrit3;
boot_crit_np = bootcrit_np1||bootcrit_np2||bootcrit_np3;
Stat_crit_values[,1] = Statistics`;
Stat_crit_values[,2] = asympt_crit`;
Stat_crit_values[,3] = boot_crit`;
Stat_crit_values [,4] = boot_crit_np`;
statrownames = { "g1" "g2" "g3" };
statcolnames = { "Statistic" "Asymptotic" "Bootstrap" "NP_Bootstrap"};
*/

*** THE MC-SIMULATION STARTS HERE;

start MCsimulation(RRgls, RRboot, RRboot_np, RRasympt, RR_PCSE, n, reps, nboots, R, T, P, sigmahat, rhohatsur, H, A, X, bparksres, bparksr, 
                                        asympt_crit_value, Xres, seed);

        gboot = j(nboots,1,0);
		p_gls = j(reps, 1,0);
        p_asym = j(reps, 1,0);
        p_boot = j(reps,1,0);
        p_boot_np = j(reps,1,0);
		p_PCSE = j(reps,1,0);
        nT = n*T;

        do k = 1 to reps;
				Testsim = 0;
				Testsimr = 0;
				Teststat = 0;

			do until (Testsim=1&Testsimr=1&Teststat=1);
                v = normal(j(n,T,seed));
                e = H`*v;
                u = j(n, T, 0);
				InvA = ginv(A);
				*run MyInv2(InvA, A);
                u[,1] = InvA*e[,1];
                do i = 2 to T;
                        u[,i] = rhohatsur*u[,i-1] + e[,i];
                end;
                ysim = X*bparksres + shape(u,nT, 1);

				*** GLS based test;
				invsigmahat = inv(sigmahat);
				InvOmega = P`*(invsigmahat@I(t))*P ;
				bsimgls = inv(X`*InvOmega*X)*(X`*InvOmega*ysim);
				vcovbsimgls = inv(X`*InvOmega*X);
				run MyInv1(Invcovbsimgls, R*vcovbsimgls*R`);
				g_simgls = (R*bsimgls)`*Invcovbsimgls*(R*bsimgls);
				if g_simgls >= asympt_crit_value then p_gls[k] = 1;




                run ParksSur(rhohatsim, Testsim, sigmahatsim, Hsim, Asim, Psim, bparkssim, vcovbparksim, T, ysim, X);
                        run MyInv1(Invbsim, R*vcovbparksim*R`);
                        g_sim = (R*bparkssim)`*Invbsim*(R*bparkssim);
                        run ParksSur(rhohatr, Testsimr, sigmahatr, Hr, Ar, Pr, bparksr, vcovbparksr, T, ysim, Xres);
                                *A_boot = Ar;
                                *H_boot = Hr;
                                *rhohat_boot = rhohatr; 
                                run Statistic(Teststat, bootstrap_crit_value, bootstrap_crit_value_np, n, nboots, R, T,
												ysim, X, bparksres, asympt_crit_value, Xres, Ar, Hr, rhohatr, bparksr, seed);
			end;
                                if g_sim > asympt_crit_value then p_asym[k] = 1;
                                if g_sim > bootstrap_crit_value then p_boot[k] = 1;
                                if g_sim > bootstrap_crit_value_np then p_boot_np[k] = 1;       

			run PCSE(vstarsim, bstarsim, covbsim, X, ysim, n, t);
			gsimPCSE = (R*bstarsim)`*inv(R*covbsim*R`)*(R*bstarsim);
			if gsimPCSE >= asympt_crit_value then p_PCSE[k] = 1;

        end;
        RRgls = p_gls[:,]; 
        RRboot = p_boot[:,];
        RRboot_np = p_boot_np[:,];
        RRasympt = p_asym[:,];
		RR_PCSE = p_PCSE[:,];

finish MCsimulation;

*run MCsimulation(RRboot1, RRboot_np1, RRasympt1, RR_PCSE1, n, reps, nboots, R, T, rhohatsur, H, A, X, bparksres, bparksr, 
                                        asympt_crit_value1, Xres, seed);
*run MCsimulation(RRboot2, RRboot_np2, RRasympt2, RR_PCSE2, n, reps, nboots, R2, T, rhohatsur, H, A, X, bparksres2, bparksr2, 
                                        asympt_crit_value2, Xres2, seed);
run MCsimulation(RRgls3, RRboot3, RRboot_np3, RRasympt3, RR_PCSE3, n, reps, nboots, R3, T, P, sigmahat, rhohatsur, H, A, X, bparksres3, bparksr3, 
                                        asympt_crit_value3, Xres3, seed);
/*
RRate = j(3,4,0);
RRate[1,] = RRboot1||RRboot_np1||RRasympt1||RR_PCSE1;
RRate[2,] = RRboot2||RRboot_np2||RRasympt2||RR_PCSE2;
RRate[3,] = RRboot3||RRboot_np3||RRasympt3||RR_PCSE3;
*/
statistics = {"g3"};
columns = {"Bootstrap" "NP_Bootstrap" "Asymptotic"};
columns2 = {"N" "T" "GLS" "Asymptotic" "Bootstrap" "NP_Bootstrap" "PCSE"};

*sderr = sqrt(vecdiag (vcovbparks));

*print bparks sderr;
*print bparks bparksres bparksres2 bparksres3;
*print rhohatsur;

*RRate2 = n||T||RRgls1||RRasympt1||RRboot1||RRboot_np1||RR_PCSE1;
*RRate2 = n||T||RRgls2||RRasympt2||RRboot2||RRboot_np2||RR_PCSE2;
RRate3 = n||T||RRgls3||RRasympt3||RRboot3||RRboot_np3||RR_PCSE3;

*print Stat_crit_values[rowname = statrownames colname = statcolnames];
print RRate3[rowname = statistics colname = columns2];

finish;
run program;
quit;

run;

